R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

source("tianfengRwrappers.R")

library(xgboost)

载入程辑包：‘xgboost’

The following object is masked from ‘package:plotly’:

    slice

The following object is masked from ‘package:dplyr’:

    slice

library(Matrix)

载入程辑包：‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

library(mclust)

    __  ___________    __  _____________
   /  |/  / ____/ /   / / / / ___/_  __/
  / /|_/ / /   / /   / / / /\__ \ / /   
 / /  / / /___/ /___/ /_/ /___/ // /    
/_/  /_/\____/_____/\____//____//_/    version 5.4.9
Type 'citation("mclust")' for citing this R package in publications.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
Registered S3 method overwritten by 'cli':
  method     from    
  print.boxx spatstat
─ Attaching packages ──────────────────────────────── tidyverse 1.3.1 ─
✓ tibble  3.1.5     ✓ stringr 1.4.0
✓ readr   2.1.2     ✓ forcats 0.5.1
✓ purrr   0.3.4     
─ Conflicts ────────────────────────────────── tidyverse_conflicts() ─
x tibble::as_data_frame()  masks igraph::as_data_frame(), dplyr::as_data_frame()
x Biobase::combine()       masks BiocGenerics::combine(), dplyr::combine()
x purrr::compose()         masks igraph::compose()
x igraph::crossing()       masks tidyr::crossing()
x Matrix::expand()         masks tidyr::expand()
x plotly::filter()         masks dplyr::filter(), stats::filter()
x widgetTools::funs()      masks dplyr::funs()
x igraph::groups()         masks plotly::groups(), dplyr::groups()
x dplyr::lag()             masks stats::lag()
x purrr::map()             masks mclust::map()
x Matrix::pack()           masks tidyr::pack()
x BiocGenerics::Position() masks ggplot2::Position(), base::Position()
x purrr::simplify()        masks igraph::simplify(), clusterProfiler::simplify()
x xgboost::slice()         masks plotly::slice(), dplyr::slice()
x Matrix::unpack()         masks tidyr::unpack()

ds2 <- readRDS("ds2.rds")
Idents(ds2) <- ds2$seurat_clusters

Error in validObject(.Object) : 
  invalid class “Seurat” object: invalid object for slot "project.name" in class "Seurat": got class "pairlist", should be or extend class "character"

数值化

ds2训练分类器

ds2_data <- get_data_table(ds2, highvar = F, type = "data")

Error in get_data_table(ds2, highvar = F, type = "data") : 
  没有"get_data_table"这个函数

importance <- xgb.importance(colnames(ds2_train), model = bst_model)
head(importance)
xgb.ggplot.importance(head(importance,20),n_clusters = 1) + theme_bw()+theme(
    axis.title.x = element_text(size = 15), axis.text.x = element_text(size = 12, colour = "black"),
    axis.title.y = element_text(size = 15), axis.text.y = element_text(size = 12, colour = "black"),
    legend.text = element_text(size = 20), legend.title = element_blank(), panel.grid = element_blank())

ds2 -> ds1

ARI = 0.1695417

Idents(ds1) <- ds1$seurat_clusters
temp <- get_data_table(ds1, highvar = F, type = "data")
ds1_data <- matrix(data=0, nrow = length(rownames(ds2_data)), ncol = length(colnames(temp)), 
                   byrow = FALSE, dimnames = list(rownames(ds2_data),colnames(temp)))
for(i in intersect(rownames(ds2_data), rownames(temp))){
  ds1_data[i,] <- temp[i,]
}
rm(temp)
ds1_label <- as.numeric(as.character(Idents(ds1)))
colnames(ds1_data) <- NULL
ds1_test_data <- list(data = t(as(ds1_data,"dgCMatrix")), label = ds1_label)
ds1_test <- xgb.DMatrix(data = ds1_test_data$data,label = ds1_test_data$label)

#预测结果
predict_ds1_test <- predict(bst_model, newdata = ds1_test)

predict_prop_ds1 <- matrix(data=predict_ds1_test, nrow = length(levels(Idents(ds2))), 
                           ncol = ncol(ds1), byrow = FALSE, 
                           dimnames = list(levels(Idents(ds2)),colnames(ds1)))

## 得到分群结果
ds1_res <- apply(predict_prop_ds1,2,func,rownames(predict_prop_ds1))
adjustedRandIndex(ds1_res, ds1_test_data$label)

Idents(ds1) <- factor(ds1_res,levels = c(0:4))
umapplot(ds1)
ds1$supclustering <- Idents(ds1) #保存监督聚类结果

数值化地投射回umap

embedding <- FetchData(object = ds1, vars = c("UMAP_1", "UMAP_2"))
embedding <- cbind(embedding, t(predict_prop_ds1))

ggobj <- ggplot() +
  geom_point(data = embedding[embedding$`0`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `0`), shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('0', low = "#FFFFFF00", high = "#6dc0a6") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`1`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `1`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('1', low = "#FFFFFF00", high = "#e2b398") +
   new_scale("color") +
    geom_point(data = embedding[embedding$`2`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `2`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('2', low = "#FFFFFF00", high = "#e2a2ca") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`3`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `3`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('3', low = "#FFFFFF00", high = "#d1eba8") +
   new_scale("color") +
      geom_point(data = embedding[embedding$`4`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `4`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('4', low = "#FFFFFF00", high = "#b1d6fb") +
    new_scale("color") +
        xlab("UMAP 1") + ylab("UMAP 2")  +
        theme(axis.line = element_line(arrow = arrow(length = unit(0.2, "cm")))) +
        scale_y_continuous(breaks = NULL) +
        scale_x_continuous(breaks = NULL) + 
  theme(panel.background = element_blank(), panel.grid = element_blank(), legend.position = "bottom")
ggsave("pre_ds1_umap.svg",device = svg,plot = ggobj,height = 10,width = 10)

ds2 -> ds0

ARI = 0.6664657

Idents(ds0) <- ds0$seurat_clusters
temp <- get_data_table(ds0, highvar = F, type = "data")
ds0_data <- matrix(data=0, nrow = length(rownames(ds2_data)), ncol = length(colnames(temp)), 
                   byrow = FALSE, dimnames = list(rownames(ds2_data),colnames(temp)))
for(i in intersect(rownames(ds2_data), rownames(temp))){
  ds0_data[i,] <- temp[i,]
}
rm(temp)
ds0_label <- as.numeric(as.character(Idents(ds0)))
colnames(ds0_data) <- NULL
ds0_test_data <- list(data = t(as(ds0_data,"dgCMatrix")), label = ds0_label)
ds0_test <- xgb.DMatrix(data = ds0_test_data$data,label = ds0_test_data$label)

#预测结果

predict_ds0_test <- predict(bst_model, newdata = ds0_test)

predict_prop_ds0 <- matrix(data=predict_ds0_test, nrow = length(levels(Idents(ds2))), 
                           ncol = ncol(ds0), byrow = FALSE, 
                           dimnames = list(levels(Idents(ds2)),colnames(ds0)))

## 得到分群结果
ds0_res <- apply(predict_prop_ds0,2,func,rownames(predict_prop_ds0))
adjustedRandIndex(ds0_res, ds0_test_data$label)
Idents(ds0) <- factor(ds0_res,levels = c(0:4))
umapplot(ds0)
ds0$supclustering <- Idents(ds0) #保存监督聚类结果

embedding <- FetchData(object = ds0, vars = c("UMAP_1", "UMAP_2"))
embedding <- cbind(embedding, t(predict_prop_ds0))

ggobj <- ggplot() +
  geom_point(data = embedding[embedding$`0`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `0`), shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('0', low = "#FFFFFF00", high = "#6dc0a6") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`1`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `1`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('1', low = "#FFFFFF00", high = "#e2b398") +
   new_scale("color") +
    geom_point(data = embedding[embedding$`2`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `2`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('2', low = "#FFFFFF00", high = "#e2a2ca") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`3`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `3`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('3', low = "#FFFFFF00", high = "#d1eba8") +
   new_scale("color") +
      geom_point(data = embedding[embedding$`4`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `4`),shape=16, size = 3, alpha=0.5) + 
  scale_color_gradient('4', low = "#FFFFFF00", high = "#b1d6fb") +
    new_scale("color") +
        xlab("UMAP 1") + ylab("UMAP 2")  +
        theme(axis.line = element_line(arrow = arrow(length = unit(0.2, "cm")))) +
        scale_y_continuous(breaks = NULL) +
        scale_x_continuous(breaks = NULL) + 
  theme(panel.background = element_blank(), panel.grid = element_blank(), legend.position = "bottom")
ggsave("pre_ds0_umap.svg",device = svg,plot = ggobj,height = 10,width = 10)

PA -> AC

Idents(ds2_PA) <- ds2_PA$seurat_clusters
selected_features <- read.csv("./datatable/selected_features.csv", stringsAsFactors = F)
selected_features <- selected_features$x
PA_data <- get_data_table(ds2_PA, highvar = F, type = "data")
PA_data <- PA_data[selected_features,]
PA_label <- as.numeric(as.character(Idents(ds2_PA)))
colnames(PA_data) <- NULL

PA_train_data <- list(data = t(as(PA_data,"dgCMatrix")), label = PA_label)
PA_train <- xgb.DMatrix(data = PA_train_data$data,label = PA_train_data$label)
xgb_param <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(ds2_PA))),
                  objective = "multi:softprob", eval_metric = 'mlogloss')

bst_model <- xgb.train(xgb_param, PA_train, nrounds = 100, verbose = 0)

Idents(ds2_AC) <- ds2_AC$seurat_clusters
AC_data <- get_data_table(ds2_AC, highvar = F, type = "data")
AC_data <- AC_data[selected_features,]
AC_label <- as.numeric(as.character(Idents(ds2_AC)))
colnames(AC_data) <- NULL
AC_test_data <- list(data = t(as(AC_data,"dgCMatrix")), label = AC_label)
AC_test <- xgb.DMatrix(data = AC_test_data$data,label = AC_test_data$label)

#预测结果
predict_prop_AC <-predict(bst_model, newdata = AC_test) %>%
 matrix(nrow = length(levels(Idents(ds2_PA))), 
                           ncol = ncol(ds2_AC), byrow = FALSE, 
                           dimnames = list(levels(Idents(ds2_PA)),colnames(ds2_AC)))
AC_res <- apply(predict_prop_AC,2,func,rownames(predict_prop_AC))

confuse_matrix1 <- table(AC_test_data$label, AC_res, dnn=c("true","pre"))
sankey_plot(confuse_matrix1,session = "PAtoAC")

Idents(ds2_AC) <- factor(AC_res,levels = c(0:2))
umapplot(ds2_AC)

embedding <- FetchData(object = ds2_AC, vars = c("UMAP_1", "UMAP_2"))
embedding <- cbind(embedding, t(predict_prop_AC))

ggobj <- ggplot() +
  geom_point(data = embedding[embedding$`0`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `0`), shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('0', low = "#FFFFFF00", high = "#6dc0a6") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`1`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `1`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('1', low = "#FFFFFF00", high = "#e2b398") +
   new_scale("color") +
    geom_point(data = embedding[embedding$`2`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `2`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('2', low = "#FFFFFF00", high = "#e2a2ca") +
        xlab("UMAP 1") + ylab("UMAP 2")  +
        theme(axis.line = element_line(arrow = arrow(length = unit(0.2, "cm")))) +
        scale_y_continuous(breaks = NULL) +
        scale_x_continuous(breaks = NULL) + 
  theme(panel.background = element_blank(), panel.grid = element_blank(), legend.position = "bottom")
ggsave("ds2_PAtoAC_umap.svg",device = svg,plot = ggobj,height = 8,width = 8)

AC to PA

Idents(ds2_AC) <- ds2_AC$seurat_clusters
selected_features <- read.csv("./datatable/selected_features.csv", stringsAsFactors = F)
selected_features <- selected_features$x
AC_data <- get_data_table(ds2_AC, highvar = F, type = "data")
AC_data <- AC_data[selected_features,]
AC_label <- as.numeric(as.character(Idents(ds2_AC)))
colnames(AC_data) <- NULL

AC_train_data <- list(data = t(as(AC_data,"dgCMatrix")), label = AC_label)
AC_train <- xgb.DMatrix(data = AC_train_data$data,label = AC_train_data$label)
xgb_ACram <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(ds2_AC))),
                  objective = "multi:softprob", eval_metric = 'mlogloss')

bst_model <- xgb.train(xgb_ACram, AC_train, nrounds = 100, verbose = 0)

Idents(ds2_PA) <- factor(ds2_PA$seurat_clusters,levels = c(0,1,2))

PA_data <- get_data_table(ds2_PA, highvar = F, type = "data")
PA_data <- PA_data[selected_features,]
PA_label <- as.numeric(as.character(Idents(ds2_PA)))
colnames(PA_data) <- NULL
PA_test_data <- list(data = t(as(PA_data,"dgCMatrix")), label = PA_label)
PA_test <- xgb.DMatrix(data = PA_test_data$data,label = PA_test_data$label)

#预测结果
predict_prop_PA <-predict(bst_model, newdata = PA_test) %>%
 matrix(nrow = length(levels(Idents(ds2_AC))), 
                           ncol = ncol(ds2_PA), byrow = FALSE, 
                           dimnames = list(levels(Idents(ds2_AC)),colnames(ds2_PA)))
PA_res <- apply(predict_prop_PA,2,func,rownames(predict_prop_PA))

confuse_matrix1 <- table(PA_test_data$label, PA_res, dnn=c("true","pre"))
sankey_plot(confuse_matrix1,session = "ACtoPA")

Idents(ds2_PA) <- factor(PA_res)
umapplot(ds2_PA)

embedding <- FetchData(object = ds2_PA, vars = c("UMAP_1", "UMAP_2"))
embedding <- cbind(embedding, t(predict_prop_PA))

ggobj <- ggplot() +
  geom_point(data = embedding[embedding$`0`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `0`), shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('0', low = "#FFFFFF00", high = "#6dc0a6") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`1`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `1`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('1', low = "#FFFFFF00", high = "#e2b398") +
   new_scale("color") +
    geom_point(data = embedding[embedding$`2`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `2`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('2', low = "#FFFFFF00", high = "#e2a2ca") +
     new_scale("color") +
    geom_point(data = embedding[embedding$`3`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `3`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('3', low = "#FFFFFF00", high = "#d1eba8") +
        xlab("UMAP 1") + ylab("UMAP 2")  +
        theme(axis.line = element_line(arrow = arrow(length = unit(0.2, "cm")))) +
        scale_y_continuous(breaks = NULL) +
        scale_x_continuous(breaks = NULL) + 
  theme(panel.background = element_blank(), panel.grid = element_blank(), legend.position = "bottom")
ggsave("ds2_ACtoPA_umap.svg",device = svg,plot = ggobj,height = 8,width = 8)

在ds0上训练

Idents(ds0) <- ds0$seurat_clusters
ds0_data <- get_data_table(ds0, highvar = F, type = "data")
ds0_label <- as.numeric(as.character(Idents(ds0)))

index <- c(1:dim(ds0_data)[2]) %>% sample(ceiling(0.3*dim(ds0_data)[2]), replace = F, prob = NULL)
colnames(ds0_data) <- NULL

ds0_train_data <- list(data = t(as(ds0_data[,-index],"dgCMatrix")), label = ds0_label[-index])
ds0_test_data <- list(data = t(as(ds0_data[,index],"dgCMatrix")), label = ds0_label[index])

ds0_train <- xgb.DMatrix(data = ds0_train_data$data,label = ds0_train_data$label)
ds0_test <- xgb.DMatrix(data = ds0_test_data$data,label = ds0_test_data$label)

watchlist <- list(train = ds0_train, eval = ds0_test)
xgb_param <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(ds0))),
                  objective = "multi:softprob", eval_metric = 'mlogloss')

bst_model <- xgb.train(xgb_param, ds0_train, nrounds = 100, watchlist, verbose = 0)
saveRDS(bst_model, "ds0_model.rds")
eval_loss <- bst_model[["evaluation_log"]][["eval_mlogloss"]]
plot_ly(data.frame(eval_loss), x = c(1:100), y = eval_loss) %>% 
  add_trace(type = "scatter", mode = "markers+lines", 
            marker = list(color = "black", line = list(color = "#1E90FFC7", width = 1)),
            line = list(color = "#1E90FF80", width = 2)) %>% 
  layout(xaxis = list(title = "epoch"),yaxis = list(title = "eval_mlogloss"))

importance <- xgb.importance(colnames(ds0_train), model = bst_model)
head(importance)
xgb.ggplot.importance(head(importance,20),n_clusters = 1) + theme_bw()+theme(
    axis.title.x = element_text(size = 15), axis.text.x = element_text(size = 8, colour = "black"),
    axis.title.y = element_text(size = 15), axis.text.y = element_text(size = 12, colour = "black"),
    legend.text = element_text(size = 20), legend.title = element_blank(), panel.grid = element_blank())
write.csv(importance, "./datatable/ds0_features.csv", row.names = F)
multi_featureplot(head(importance,9)$Feature, ds0, labels = "") 

ds0 -> ds2

ARI = 0.4015002

Idents(ds2) <- ds2$seurat_clusters 
temp <- get_data_table(ds2, highvar = F, type = "data")
ds2_data <- matrix(data=0, nrow = length(rownames(ds0_data)), ncol = length(colnames(temp)), 
                   byrow = FALSE, dimnames = list(rownames(ds0_data),colnames(temp)))
for(i in intersect(rownames(ds2_data), rownames(temp))){
  ds2_data[i,] <- temp[i,]
}
rm(temp)
ds2_label <- as.numeric(as.character(Idents(ds2)))
colnames(ds2_data) <- NULL
ds2_test_data <- list(data = t(as(ds2_data,"dgCMatrix")), label = ds2_label)
ds2_test <- xgb.DMatrix(data = ds2_test_data$data,label = ds2_test_data$label)

#预测结果

predict_ds2_test <- predict(bst_model, newdata = ds2_test)

predict_prop_ds2 <- matrix(data=predict_ds2_test, nrow = bst_model[["params"]][["num_class"]], 
                           ncol = ncol(ds2), byrow = FALSE, 
                           dimnames = list(c(0:(bst_model[["params"]][["num_class"]]-1)),colnames(ds2)))

## 得到分群结果
ds2_res <- apply(predict_prop_ds2,2,func,rownames(predict_prop_ds2))

adjustedRandIndex(ds2_res, ds2_test_data$label)
confuse_matrix1 <- table(ds2_test_data$label, ds2_res, dnn=c("true","pre"))

sankey_plot(confuse_matrix1,0:5,0:4,session = "ds0tods2")

Idents(ds2) <- factor(ds2_res,levels = c(0:5))
umapplot(ds2)

embedding <- FetchData(object = ds2, vars = c("UMAP_1", "UMAP_2"))
embedding <- cbind(embedding, t(predict_prop_ds2))

ggobj <- ggplot() +
  geom_point(data = embedding[embedding$`0`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `0`), shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('0', low = "#FFFFFF00", high = "#6dc0a6") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`1`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `1`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('1', low = "#FFFFFF00", high = "#e2b398") +
   new_scale("color") +
    geom_point(data = embedding[embedding$`2`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `2`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('2', low = "#FFFFFF00", high = "#e2a2ca") +
     new_scale("color") +
    geom_point(data = embedding[embedding$`3`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `3`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('3', low = "#FFFFFF00", high = "#d1eba8") +
     new_scale("color") +
    geom_point(data = embedding[embedding$`4`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `4`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('4', low = "#FFFFFF00", high = "#b1d6fb") +
     new_scale("color") +
    geom_point(data = embedding[embedding$`5`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `5`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('5', low = "#FFFFFF00", high = "#fd9999") +
        xlab("UMAP 1") + ylab("UMAP 2")  +
        theme(axis.line = element_line(arrow = arrow(length = unit(0.2, "cm")))) +
        scale_y_continuous(breaks = NULL) +
        scale_x_continuous(breaks = NULL) + 
  theme(panel.background = element_blank(), panel.grid = element_blank(), legend.position = "bottom")
ggsave("ds0tods2umap.svg",device = svg,plot = ggobj,height = 8,width = 8)

ds0 -> ds1

ARI = 0.2524238

Idents(ds1) <- ds1$seurat_clusters
temp <- get_data_table(ds1, highvar = F, type = "data")
ds1_data <- matrix(data=0, nrow = length(rownames(ds0_data)), ncol = length(colnames(temp)), 
                   byrow = FALSE, dimnames = list(rownames(ds0_data),colnames(temp)))
for(i in intersect(rownames(ds1_data), rownames(temp))){
  ds1_data[i,] <- temp[i,]
}
rm(temp)
ds1_label <- as.numeric(as.character(Idents(ds1)))
colnames(ds1_data) <- NULL
ds1_test_data <- list(data = t(as(ds1_data,"dgCMatrix")), label = ds1_label)
ds1_test <- xgb.DMatrix(data = ds1_test_data$data,label = ds1_test_data$label)

#预测结果

predict_ds1_test <- predict(bst_model, newdata = ds1_test)

predict_prop_ds1 <- matrix(data=predict_ds1_test, nrow = bst_model[["params"]][["num_class"]], 
                           ncol = ncol(ds1), byrow = FALSE, 
                           dimnames = list(c(0:(bst_model[["params"]][["num_class"]]-1)),colnames(ds1)))

## 得到分群结果
ds1_res <- apply(predict_prop_ds1,2,func,rownames(predict_prop_ds1))
adjustedRandIndex(ds1_test_data$label, ds1_res)
Idents(ds1) <- factor(ds1_res,levels = c(0:5))
umapplot(ds1)
# umapplot(ds1,group.by = "Classification1")
confuse_matrix <- table(ds1_test_data$label, ds1_res, dnn=c("true","pre"))
sankey_plot(confuse_matrix,c(0:4),c(0:4),session = "ds0tods1")

embedding <- FetchData(object = ds1, vars = c("UMAP_1", "UMAP_2"))
embedding <- cbind(embedding, t(predict_prop_ds1))

ggobj <- ggplot() +
  geom_point(data = embedding[embedding$`0`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `0`), shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('0', low = "#FFFFFF00", high = "#6dc0a6") +
  new_scale("color") +
    geom_point(data = embedding[embedding$`1`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `1`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('1', low = "#FFFFFF00", high = "#e2b398") +
   new_scale("color") +
    geom_point(data = embedding[embedding$`2`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `2`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('2', low = "#FFFFFF00", high = "#e2a2ca") +
     new_scale("color") +
    geom_point(data = embedding[embedding$`3`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `3`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('3', low = "#FFFFFF00", high = "#d1eba8") +
     new_scale("color") +
    geom_point(data = embedding[embedding$`4`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `4`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('4', low = "#FFFFFF00", high = "#b1d6fb") +
     new_scale("color") +
    geom_point(data = embedding[embedding$`5`>0.1,], 
             aes(x = UMAP_1, y = UMAP_2, color = `5`),shape=16, size = 2, alpha=0.5) + 
  scale_color_gradient('5', low = "#FFFFFF00", high = "#fd9999") +
        xlab("UMAP 1") + ylab("UMAP 2")  +
        theme(axis.line = element_line(arrow = arrow(length = unit(0.2, "cm")))) +
        scale_y_continuous(breaks = NULL) +
        scale_x_continuous(breaks = NULL) + 
  theme(panel.background = element_blank(), panel.grid = element_blank(), legend.position = "bottom")
ggsave("ds0tods1umap.svg",device = svg,plot = ggobj,height = 8,width = 8)

ARI 和聚类数的关系

lym_ds2 <- readRDS("./lym_ds2.rds")

Idents(lym_ds2) <- lym_ds2$conditions
lym_ds2_AC <- subset(lym_ds2, idents = "AC")
lym_ds2_PA <- subset(lym_ds2, idents = "PA")

set.seed(7)
lym_PA_data <- get_data_table(lym_ds2_PA, highvar = T, type = "data")
lym_AC_data <- get_data_table(lym_ds2_AC, highvar = T, type = "data")
res <- list()

for(reso in seq(0.05,0.3,0.05))
{
  lym_ds2_AC <- lym_ds2_AC %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = reso)
  lym_ds2_PA <- lym_ds2_PA %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = reso)

  lym_PA_label <- as.numeric(as.character(Idents(lym_ds2_PA)))

  index <- c(1:dim(lym_PA_data)[2]) %>% sample(ceiling(0.3*dim(lym_PA_data)[2]), replace = F, prob = NULL)
  colnames(lym_PA_data) <- NULL
  lym_PA_train_data <- list(data = t(as(lym_PA_data[,-index],"dgCMatrix")), label = lym_PA_label[-index])
  lym_PA_test_data <- list(data = t(as(lym_PA_data[,index],"dgCMatrix")), label = lym_PA_label[index])
  
  lym_PA_train <- xgb.DMatrix(data = lym_PA_train_data$data,label = lym_PA_train_data$label)
  lym_PA_test <- xgb.DMatrix(data = lym_PA_test_data$data,label = lym_PA_test_data$label)
  
  watchlist <- list(train = lym_PA_train, eval = lym_PA_test)
  xgb_param <- list(eta = 0.2, max_depth = 6, 
                    subsample = 0.6,  num_class = length(table(Idents(lym_ds2_PA))),
                    objective = "multi:softmax", eval_metric = 'mlogloss')
  
  bst_model <- xgb.train(xgb_param, lym_PA_train, nrounds = 100, watchlist, verbose = 0)
  
  
  lym_AC_label <- as.numeric(as.character(Idents(lym_ds2_AC)))
  colnames(lym_AC_data) <- NULL
  lym_AC_test_data <- list(data = t(as(lym_AC_data,"dgCMatrix")), label = lym_AC_label)
  lym_AC_test <- xgb.DMatrix(data = lym_AC_test_data$data,label = lym_AC_test_data$label)
  predict_lym_AC_test <- round(predict(bst_model, newdata = lym_AC_test))
  
  res <- append(x = res,values = adjustedRandIndex(predict_lym_AC_test, lym_AC_test_data$label))
}
res

ARI 和聚类数的关系

 ds2 <- readRDS("./ds2.rds")

Idents( ds2) <-  ds2$conditions
 ds2_AC <- subset( ds2, idents = "AC")
 ds2_PA <- subset( ds2, idents = "PA")

set.seed(7)
 PA_data <- get_data_table( ds2_PA, highvar = T, type = "data")
 AC_data <- get_data_table( ds2_AC, highvar = T, type = "data")
ds2_res <- list()

for(reso in seq(0.05,0.3,0.05))
{
   ds2_AC <-  ds2_AC %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = reso)
   ds2_PA <-  ds2_PA %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = reso)

   PA_label <- as.numeric(as.character(Idents( ds2_PA)))

  index <- c(1:dim( PA_data)[2]) %>% sample(ceiling(0.3*dim( PA_data)[2]), replace = F, prob = NULL)
  colnames( PA_data) <- NULL
   PA_train_data <- list(data = t(as( PA_data[,-index],"dgCMatrix")), label =  PA_label[-index])
   PA_test_data <- list(data = t(as( PA_data[,index],"dgCMatrix")), label =  PA_label[index])
  
   PA_train <- xgb.DMatrix(data =  PA_train_data$data,label =  PA_train_data$label)
   PA_test <- xgb.DMatrix(data =  PA_test_data$data,label =  PA_test_data$label)
  
  watchlist <- list(train =  PA_train, eval =  PA_test)
  xgb_param <- list(eta = 0.2, max_depth = 6, 
                    subsample = 0.6,  num_class = length(table(Idents( ds2_PA))),
                    objective = "multi:softmax", eval_metric = 'mlogloss')
  
  bst_model <- xgb.train(xgb_param,  PA_train, nrounds = 100, watchlist, verbose = 0)
  
  
   AC_label <- as.numeric(as.character(Idents( ds2_AC)))
  colnames( AC_data) <- NULL
   AC_test_data <- list(data = t(as( AC_data,"dgCMatrix")), label =  AC_label)
   AC_test <- xgb.DMatrix(data =  AC_test_data$data,label =  AC_test_data$label)
  predict_AC_test <- round(predict(bst_model, newdata = AC_test))
  
  ds2_res <- append(x = ds2_res,values = adjustedRandIndex(predict_AC_test,  AC_test_data$label))
}
ds2_res

绘图

plot_ly(data, x = ~resolution) %>% 
  add_trace(y = ~lym_ARI, name = 'lym', mode = 'lines+markers',size = 4, markers = list(color = "#158aff"),
            line = list(color = '#b1d6fb', width = 4), type = 'scatter') %>%
  add_trace(y = ~SMC_ARI, name = 'SMC', mode = 'lines+markers',size = 4, markers = list(color = "#ff2121"), 
            line = list(color = '#e2a2ca', width = 4), type = 'scatter')%>%
  layout(xaxis = list(zerolinecolor = '#ffffff', zerolinewidth = 2, range = c(0, 0.32),
                       gridcolor = 'ffff', dtick=0.05, title = "resolution"),
         yaxis = list(zerolinecolor = '#ffffff', zerolinewidth = 2, range = c(0.2, 0.9),
                       gridcolor = 'ffff', dtick=0.1,  title = "ARI")) %>%
  layout(font = list(family = "Arial", size = 20, color = "black"))

Warning: 'scatter' objects don't have these attributes: 'markers'
Valid attributes include:
'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isG [... truncated]
Warning: 'scatter' objects don't have these attributes: 'markers'
Valid attributes include:
'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isG [... truncated]
Warning: 'scatter' objects don't have these attributes: 'markers'
Valid attributes include:
'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isG [... truncated]
Warning: 'scatter' objects don't have these attributes: 'markers'
Valid attributes include:
'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isG [... truncated]

解释性

library(DALEX)

Welcome to DALEX (version: 2.4.1).
Find examples and detailed introduction at: http://ema.drwhy.ai/

library(modelStudio)
source("tianfengRwrappers.R")

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
载入需要的程辑包：dplyr

载入程辑包：‘dplyr’

The following object is masked from ‘package:DALEX’:

    explain

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

载入需要的程辑包：reticulate
载入需要的程辑包：tidyr

载入程辑包：‘MySeuratWrappers’

The following objects are masked from ‘package:Seurat’:

    DimPlot, DoHeatmap, LabelClusters, RidgePlot, VlnPlot


载入程辑包：‘cowplot’

The following object is masked from ‘package:ggpubr’:

    get_legend

载入需要的程辑包：viridisLite

载入程辑包：‘reshape2’

The following object is masked from ‘package:tidyr’:

    smiths

NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
      Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
      if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow

Registered S3 method overwritten by 'enrichplot':
  method               from
  fortify.enrichResult DOSE
clusterProfiler v3.14.3  For help: https://guangchuangyu.github.io/software/clusterProfiler

If you use clusterProfiler in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.
Registering fonts with R

载入程辑包：‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

载入需要的程辑包：Biobase
载入需要的程辑包：BiocGenerics
载入需要的程辑包：parallel

载入程辑包：‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply,
    parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated,
    eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match,
    mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames,
    sapply, setdiff, sort, table, tapply, union, unique, unsplit, which, which.max, which.min

Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

载入需要的程辑包：e1071

载入程辑包：‘widgetTools’

The following object is masked from ‘package:dplyr’:

    funs


载入程辑包：‘DynDoc’

The following object is masked from ‘package:BiocGenerics’:

    path


载入程辑包：‘DT’

The following object is masked from ‘package:Seurat’:

    JS

========================================
circlize version 0.4.13
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
  in R. Bioinformatics 2014.

This message can be suppressed by:
  suppressPackageStartupMessages(library(circlize))
========================================

载入需要的程辑包：grid
========================================
ComplexHeatmap version 2.2.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.
========================================


载入程辑包：‘ComplexHeatmap’

The following object is masked from ‘package:plotly’:

    add_heatmap

ds2 <- readRDS("ds2.rds")

# 建立解释器
explainer_xgb <- DALEX::explain(bst_model, data = ds2_train_data$data[1:200,], y = ds2_train_data$label[1:200], label = "XGBoost")

# modelStudio::modelStudio(explainer = explainer_xgb, new_observation = ds2_train_data$data[1:2,] )
p1 <- DALEX::variable_profile(explainer_xgb, variable = "ACKR1", type = "accumulated")
plot(p1)
p2 <- single_variable(explainer_xgb, rownames(ds2_data)[1:10], type = "pdp")
plot(p2)
# p3 <- variable_importance(explainer_xgb, variable =rownames(ds2_data)[1:10],loss_function = loss_root_mean_square)

reduced_ds2model <- readRDS("reduced_ds2model.rds")
p1 <- xgb.plot.tree(feature_names = rownames(ds2_data), model = reduced_ds2model,render = F)
DiagrammeR::render_graph(p1)
saveRDS(p1,"ds2_tree.rds")

p2 <- xgb.ggplot.shap.summary(data = ds2_train_data$data, features = rownames(ds2_data)[1:10], model = reduced_ds2model)

xgb.ggplot.deepness(model = reduced_ds2model, which = c("2x1", "max.depth", "med.depth", "med.weight"))

shapvalues <- xgb.plot.shap(data = ds2_train_data$data, features = rownames(ds2_data)[1:10], model = reduced_ds2model,plot = F)

library(SHAPforxgboost)

shap_contrib <- predict(reduced_ds2model, ds2_train_data$data, predcontrib = TRUE)
shap_contrib <- data.frame(shap_contrib[[2]]) ##提取fbm对应的分类
BIAS0 <- shap_contrib[, ncol(shap_contrib)][1]
shap_contrib[,"BIAS"] <- NULL
imp <- colMeans(abs(shap_contrib))
mean_shap_score <- imp[order(imp, decreasing = T)]

shap_values<- list(shap_score = shap_contrib, mean_shap_score = mean_shap_score, 
    BIAS0 = BIAS0)

# The ranked features by mean |SHAP|
shap_values$mean_shap_score

        MYL9          LUM          CFH         DSTN          MGP         VCAN        ACTA2          BGN       COL3A1        TIMP1       COL1A2       CRTAC1         TPM2 
0.6204572985 0.5820071027 0.4785667498 0.3556163464 0.3130181710 0.2913836150 0.2837963942 0.2679459393 0.1930709253 0.1873566211 0.1505724998 0.1290744804 0.1201026216 
       TAGLN     PPP1R14A         RGS5        RAMP1         CCL2        TIMP3       IGFBP6         PI16      COL15A1          DCN       IGFBP5          AGT         APOE 
0.1003337457 0.0857010119 0.0852767693 0.0795046269 0.0757562707 0.0701276488 0.0678183694 0.0649776709 0.0631916199 0.0611326082 0.0601633480 0.0580075152 0.0578228089 
       MYH11       CHRDL1          A2M         CTSZ        SPARC         CST3        NUPR1        DAAM1       MT.CO1     SERPINE1          GSN      SLC25A4         TCF4 
0.0541330588 0.0532325922 0.0530933256 0.0513724617 0.0481444628 0.0478056960 0.0443659350 0.0438817839 0.0418114645 0.0413728777 0.0404521745 0.0403358532 0.0403044761 
       CALD1          CFD        ABCC9        ADIRF        SULF1         CTGF        MFAP4      HCFC1R1         TPI1       TMSB10         CTSC        HSPB1        RPS3A 
0.0401368241 0.0387895560 0.0384200689 0.0365435653 0.0353001928 0.0341662372 0.0341286377 0.0332152959 0.0331589190 0.0328425296 0.0325752880 0.0301465333 0.0298549392 
      TUBA1A           C7      SELENOW        CCL19      RARRES2        HSPA5       TMSB4X         ACTB         NPC2        CSRP1       LHFPL6         TGM2        HLA.B 
0.0297061329 0.0277430090 0.0276101055 0.0271826681 0.0269783002 0.0259206774 0.0253589819 0.0253233452 0.0253218930 0.0252741448 0.0244636626 0.0241604054 0.0228003503 
       MYH10     C11orf96        PLAC9      TINAGL1        ZFP36        AEBP1        TXNIP         CTSK      SPARCL1        FBLN1          OMD         FLNA        EPAS1 
0.0227170044 0.0225326449 0.0223703709 0.0222978286 0.0220347147 0.0220116509 0.0214403874 0.0209068173 0.0209057632 0.0203254536 0.0198206404 0.0196002415 0.0189674272 
        FRZB         SRGN        RPL10      TSC22D1        MEDAG       EFEMP1       EEF1A1        ITM2B       COL6A1        TFPI2          LOX       ANTXR1        RPL32 
0.0187382099 0.0185984747 0.0180381953 0.0179055931 0.0178611410 0.0178069619 0.0176033985 0.0174932844 0.0174558985 0.0173056289 0.0172437097 0.0172309449 0.0172202093 
      COL6A3      TNFSF10          CLU         MT2A        MGST1         TPT1         CD63       MT.CO2        LTBP2       IGFBP7         RGS2        MXRA7        POSTN 
0.0172045728 0.0172023475 0.0168330951 0.0165347722 0.0161459466 0.0160803324 0.0160172932 0.0159915421 0.0158408501 0.0158366315 0.0151780958 0.0149119464 0.0146892244 
        SOST        CEBPB       CCDC80       MT.ND2          CPE        PHGDH          B2M       IGFBP2         NNMT        PRKG1         GAS6        SUGCT       RPS15A 
0.0144515973 0.0144285185 0.0144071938 0.0143921291 0.0143113251 0.0141186740 0.0140907003 0.0137512746 0.0137333326 0.0137030245 0.0135413042 0.0135100512 0.0133088363 
        RPS2         APOD        DUSP1        MFAP5       OLFML3       S100A4         TPM1      TGFB1I1        EDIL3        CYR61        INHBA         MCAM      S100A10 
0.0132484229 0.0129346542 0.0126828317 0.0126363272 0.0125613754 0.0122139623 0.0120825502 0.0120522044 0.0120383867 0.0119389259 0.0119174619 0.0118887791 0.0118223191 
        BTG2        MAP1B          APP        PDE1A      TSC22D3       CXCL12         IER2       MT.ND4        RPS14        RPL34     HIST1H4C         NEXN       NFKBIZ 
0.0117749544 0.0117481146 0.0117057740 0.0115826632 0.0115442764 0.0113382007 0.0112637513 0.0112401562 0.0112064757 0.0110264808 0.0107459507 0.0106657271 0.0103559719 
      MT.ND3       LGALS1         FGF7          CD9     HLA.DPB1      S100A11        ACTN1          JUN         KLF2        PDE3A       HSPA1A        CEBPD        GNG11 
0.0101976992 0.0101770998 0.0101671011 0.0099403897 0.0099097975 0.0097435482 0.0097155518 0.0096896182 0.0096807902 0.0096555487 0.0096069540 0.0095722914 0.0095562394 
      S100A6      HLA.DRA         JAG1         IL33       TCEAL4        ADH1B          FAP          LBH         MT1E        NR2F2    TNFRSF12A        SOCS3         SAT1 
0.0095482029 0.0094487949 0.0094098613 0.0093210160 0.0092902735 0.0092513737 0.0091262500 0.0091033193 0.0090539192 0.0088978601 0.0087584394 0.0087017552 0.0086825079 
       FXYD1        CRIM1          NOV       ARID5B         IGKC         PKIG         CAV2         RPS8          FOS     NDUFA4L2         JUND       PCOLCE       TM4SF1 
0.0086061487 0.0085865603 0.0085773613 0.0083657040 0.0082623845 0.0082452222 0.0082089682 0.0079698077 0.0078723466 0.0077612101 0.0075944786 0.0075450245 0.0075007141 
       BICC1         CD59         PTMS        ACTC1         FTH1        NR4A1     HLA.DRB5         SDC2          OGN       LGALS3         PAWR      ZFP36L2          DPT 
0.0074931883 0.0074754157 0.0074396443 0.0073509889 0.0073372422 0.0071350127 0.0071314588 0.0070796633 0.0070266357 0.0070074144 0.0070027578 0.0069374098 0.0068969687 
       RPLP1          C1R     SERPINE2        CRIP2        WWTR1       COL5A2         THY1     HSP90AA1       COX7A1         ENAH       AKR1C2        CRYAB         PLS3 
0.0068921354 0.0068781974 0.0068765462 0.0068477642 0.0067674040 0.0066609383 0.0066312031 0.0065710922 0.0065462989 0.0064975395 0.0064901181 0.0064469821 0.0064368523 
       PALLD         MT1M       PIK3R1       SYNPO2       STEAP1         TYMP         SSR4      SELENOK          VIM         IER3        PLCG2        NEAT1         MT1X 
0.0063227409 0.0061759413 0.0061701969 0.0061537421 0.0061444942 0.0060946192 0.0060678631 0.0060592792 0.0060568250 0.0060264915 0.0059864846 0.0059044698 0.0058888294 
       RPS19         FHL1       IFITM1         ENO1        RCAN2          TNC       SEMA3C          LPP      ARHGDIB         BTG1        EFHD1     SH3BGRL3         LMNA 
0.0058212521 0.0057120912 0.0057067324 0.0056832581 0.0056795813 0.0056789069 0.0056038259 0.0055865127 0.0055858772 0.0055846104 0.0055603855 0.0055398404 0.0055082398 
        FHL5          FST        TSHZ2           C3       COL8A1       FKBP1A        CYTOR        VCAM1         MEG3        IER5L         BST2         TRAC        SEPT1 
0.0054868822 0.0054634564 0.0054093946 0.0053586680 0.0053346375 0.0053287444 0.0053075022 0.0052665941 0.0052164661 0.0051720644 0.0051266265 0.0051175765 0.0051172153 
     MT.ATP6     LGALS3BP          MDK        FCGRT         SSPN      SELENOM        ANXA1         NFIA        RPS18      ZFP36L1       MT.CO3       MT.CYB        SFRP1 
0.0050294289 0.0049997732 0.0049689475 0.0049372631 0.0049265895 0.0048792597 0.0048649981 0.0047039796 0.0046924844 0.0046087965 0.0045480155 0.0044777447 0.0044434913 
       LAMB1       TUBA1B         DKK3         FOSB      GADD45A        ABCA6         SMC4        ISG15       PDGFRL        HTRA1        LAMB2          ID2       IFITM3 
0.0044056135 0.0043820581 0.0042295295 0.0042172200 0.0041783916 0.0041095395 0.0040696254 0.0039991583 0.0038915954 0.0037544557 0.0037276164 0.0037212183 0.0037167951 
       MRVI1         DLX5        MXRA8        ITGA8        SFRP2         TFPI        PLPP1        RPL28        HSPH1        CCL21         TOB1       TBC1D4       COL6A2 
0.0036786448 0.0036732905 0.0036597685 0.0036233183 0.0035865804 0.0035087005 0.0034705618 0.0034511603 0.0034408128 0.0034193244 0.0033870992 0.0033180242 0.0033105706 
       FSTL1        IGLC2      COL12A1        SMOC2         STK4         ISLR        CD151      PCOLCE2        ANXA2         SOD3         LDB2          PLN      C2orf40 
0.0033015458 0.0032463733 0.0032342356 0.0031852721 0.0031582437 0.0031431667 0.0030488996 0.0030454291 0.0030252269 0.0029658435 0.0029367448 0.0029089800 0.0029067774 
       IL6ST          AXL     HLA.DRB1         CD74         SOX4       LRRC17        RPS29         DIO2      COLEC11        ZC3H6         CTSF        SLIT3         PGK1 
0.0028953885 0.0028844848 0.0028620202 0.0028240649 0.0028113534 0.0027445431 0.0026793188 0.0026665496 0.0026475506 0.0026459482 0.0026235411 0.0026144796 0.0025999836 
         ELN          FN1         FGL2        LIMA1       MICAL2          ID3      S100A16         ASPN       PDLIM3        RGS16         GPX3         CCL5        RPS12 
0.0025652139 0.0025461112 0.0025320984 0.0025270031 0.0025036490 0.0024534509 0.0024081292 0.0024048949 0.0023829600 0.0023667976 0.0023197245 0.0022136439 0.0022100411 
       SPRY1         CAV1      GUCY1B1        FABP4       LIMCH1         SNCG       IGSF10         PDK4         TBX2       NFKBIA       CALCRL         JUNB      HERPUD1 
0.0021267428 0.0021091479 0.0020465965 0.0020464108 0.0020239332 0.0019699449 0.0019670940 0.0019480049 0.0019417005 0.0018255686 0.0017568510 0.0017466342 0.0016701163 
        CREM       IGFBP3       EPSTI1         RHOB     SERPINH1         KRT8         FBN1       ANGPT2       SUCNR1        LMOD1     CRISPLD1      COL14A1      S100A13 
0.0016652647 0.0016409781 0.0016070931 0.0015988775 0.0015970256 0.0015964362 0.0014773414 0.0014578561 0.0014339205 0.0013835100 0.0013738182 0.0013687167 0.0013595937 
   MTRNR2L12       MARCKS     SERPINB1        RBPMS         LRP1     SERPING1        ARL4D        PDE5A         XAF1      ILF3.DT       C9orf3        H2AFZ      PLEKHF1 
0.0013460439 0.0012984280 0.0012954635 0.0012833999 0.0012696784 0.0012370948 0.0012229453 0.0012147770 0.0011684844 0.0011556226 0.0011503195 0.0011341533 0.0011307312 
       BIRC3       EFEMP2          C1S        VEGFC        PROCR        SRSF7         GBP1           CP       PDGFRB         AQP1          GEM        DEPP1         CD37 
0.0010614697 0.0010591948 0.0010119252 0.0009980423 0.0009882408 0.0009596119 0.0008773206 0.0008418845 0.0008366114 0.0006864486 0.0006002834 0.0005437652 0.0005223113 
   TNFRSF11B       ABI3BP         SPP1        ITLN1         CCL4       TPSAB1       S100A8       CCL4L2        TPSB2       S100A9          FTL         CPA3       JCHAIN 
0.0004900735 0.0003462739 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        GNLY        IGHA1      PLA2G2A        APOC1        CCL18        PTGDS         AREG        FABP5       CXCL10        IGHG1         SLPI         IBSP       CXCL14 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        C1QB         CCL3        IGHG3        SFRP4         MT1G        ACKR1         EDN1         C1QA          LYZ        IGHG2         MMP9        IGLC3         XCL1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        MZB1        IGHG4         XCL2         IL1B         C1QC         MT1H         IGHM        CXCL9        CCL13       CCL3L1        CCL20         CCL8     HLA.DPA1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        G0S2        CXCL3         PRG4         MMP7       RNASE1         CCL7        CXCL8        CXCL1      S100A12         HPGD      IGKV1.5       COL1A1        MS4A2 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        CTSG        CH25H         CSTB     IGHV3.49         NKG7     IGLV3.21        CCL17      SELENOP         IFNG        IL1RN         LY6H       CHI3L1        CXCL2 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       IGLL5     HLA.DQA1       CXCL13       SPINK1         FCN1     DNASE1L3       FCER1A        IGHA2      ADCYAP1         GZMB        WISP2       C2CD4B          PTN 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       TCL1A       CXCL11         TAC1        IFI27         SAA1         GZMK     SERPINB2      IGLV3.1         CD36     IGHV3.74     C15orf48        MMP12       FGFBP2 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
    SERPINF1        FCGBP         CD69         IL7R        LYVE1         CD14        RERGL        NRXN1        IGHGP         IL13          HDC         PTX3        HMOX1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        PLTP        F13A1         FBP1       TYROBP        CXCL5     APOBEC3A        S100B        IGLC7        FOLR2        CLDN5      RARRES1        SFRP5        DERL3 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
    IGKV3.20      IGKV4.1        ACKR3         STC1         CMA1     IGLV6.57          VWF         CTSL         AIF1         SELE         EREG          LTB     IGHV3.23 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       CD79A         CTSD     IGLV1.40        RAMP2     IGKV3.15        GATA2        KRT17      IGHV3.7        CYTL1         GZMA     HLA.DQA2         COMP       IL1RL1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       LYPD2         GZMH      TNFRSF4     IGLV2.14          DES        MARCO     IGHV4.39         MYOC          NTS          ID1     IGHV3.33        SOX18         RHEX 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       KRT14       ANKRD1        DUSP2        MMP13     DLX6.AS1     IGLV1.51        TRBC1        HPGDS         HHIP         IL32        KRT18       FCER1G        MS4A1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        RGCC         CD1C       STEAP4     IGHV1.24         CSF2         PPBP        PRSS2         PLAT           HP      ANGPTL7         IDO1     HIST1H1B         RBP7 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       HSPA6          KIT      IGLV2.8        KLRB1          CA4        RAMP3     SLC9A3R2        ACTG2         CTSB        CENPF        PLAUR     IGLV2.23          IL6 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        ACP5     IGKV3.11          TNF       MS4A6A        THBS1         IL10        KLRD1          SDS        IFIT2         MRC1        TRBC2     IGHV1.18         TRDC 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       CCL22       CLEC9A        FBLN2       IGFBP4        PLIN2      SLC18A2        KRT81        CXCR4        MKI67         CTSW          ID4         FMO2     HLA.DQB1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       MMP11     IGHV3.15          GAL        CCL26        PLVAP        HYAL2      CLEC10A       CRABP1         PLEK          SCT     IGKV1.39        KRT16        GPHA2 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
      FCGR3A        UBE2C         CPVL        KLRC1         IGHD         MMP1         TNXB     IGHV3.21       GPR183         IGF1    IGHV1.69D       CAVIN2         CTSS 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
    IGLV2.11         KRT7         ASPM          SLN      IGHV1.3         LGMN         MNDA         PMCH         GAS1        TIMD4        TOP2A         CST7        GDF15 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       FFAR3         CCR7         VMO1        STMN2        IFIT1        IL17A      CLEC14A      ADAMTS1         RGS1        KRT86        IGLC6     TNFRSF18       HSPA1B 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       CCL23         CSF3       CHI3L2         NEFM          CD2       LGALS2          GRP         RETN       MS4A4A         CD68     IGHV1.46         TCIM         CD5L 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
    IGLV1.44        RGS13        PTGS2         SAA2     IGKV1D.8        GREM1       COL4A1      FILIP1L        IGSF6      PLA2G2D       DNAAF1       AKAP12         TDO2 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        PRF1         CD8A   AC233755.2        GREM2     IGHV3.43      GADD45B        MMRN1         CNN1          BMX         TCN1       CYP1B1       NCCRP1       SCARA5 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       ADH1C        SPON2       NUSAP1        MS4A7        CRTAM         MMP2         CD8B     IGHV5.51        GAP43          IL2      TNFAIP6     IGLV7.43        CLIC3 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
      LILRA4         ESM1   AC020656.1        SSUH2         DKK2     IGHV4.59         PCK1          CD7       COX4I2         GJA5         SGK1     IGLV1.47     IGKV1.27 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       OBP2A        RASD1        EGFL6       INSIG1         HES1        PRELP        CSRP2      SCGB3A2         VWA1         NOX5       SPRR2F         EGR1        SOX17 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
   IGKV2D.28         EMCN         SCGN         NPPC        STMN1         CD52        C5AR1       PRSS23          LTA         HAMP        CCL11       HIGD1B     TMEM176B 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        PRND     IGHV3.11         AQP2        IFIT3     IGLV3.19        HSPG2         GJA4         TYMS           F3        KRT25        TKTL1        CRLF1       CTHRC1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       TRGC2         NEFL        HCAR3        IL1R2        CLDN1         CD1E       PECAM1        MPZL2         MAFB     C12orf75        JSRP1         HOPX        CD160 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       CD177        KRT19        DDIT4         PRPH      IGHV2.5        TREM2         MT1A        FOLR3         ACY3         PENK         RBP4        ECSCR          F10 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       MFGE8     IGHV3.48        EGFL7         LEPR        PODXL       FCGR1A         TFF3          HBB         APLN   AC233755.1         HEG1        LTBP1         CD83 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
      VPREB3         DKK1   AC092484.1        BANK1      HLA.DMA          PGF      SULT1E1      C1QTNF3          LIF   AC084871.2         CD3D       BCL2A1     CCDC102B 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       PTHLH        ABCA8      IGKV1.9   AL031316.1       CYP4B1     IGKV1.17       FCGR2A        PALMD        CARMN        KLRF1       FAM30A         MSR1        GPNMB 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
    SERPINA1     IGHV3.20       GIMAP7       ZNF331        MMP10        RSPO2     IGHV4.34         RRAD        IGLL1        HCAR2         IGHE        CALB2        CRHBP 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       SCRG1         SCG3      GUCY1A1         CYBB        EDNRB        VSIG4        HTRA3        PCLAF    LINC01781         IRF1        PRRX1        CD163        CXCR6 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
    IGHV3.30     PGM5.AS1       CLEC4C       SORBS2        KLRC2       PAPPA2         PI15     IGLV4.69          IL3         OLR1         FCN3        CCND1       CYP1A1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        SPIC    IGKV1D.39         HBA1         ADM5   IGHV1.69.2         MSMP         CD3G         MYLK        THBS2        CTLA4     ADAMDEC1       ITGBL1        TIGIT 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        LRG1      SLC40A1        CHIT1        LTC4S         CPB1       TNFSF9        GAPDH        RTKN2        CCL14          HBD       PDLIM1        IL36B         SPIB 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
      FCGR2B        NAMPT         ATF3       TSPAN7      TNFAIP3         KLF6      KIR2DL4         AGR2       ADGRL4        TIMP4       PLA2G7        OTUD1       DNAJB1 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
         LPL        CD247        GPR34         PAEP         GLUL        S100P        MTUS1          TTN         NDNF         RRM2        CDC20        BASP1       SLAMF9 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
       PTGIS        TPPP3        ITM2A        ZNF80      FAM180B        LAIR2        LTBP4        NPDC1        BIRC5        IGFL2        CD248         SCG2      IGHV4.4 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        CD70        PROK2         LST1        TPSD1     IGHV4.31         KLK1        PTTG1         SOD2        PLAC8       ZNF683        FABP3         SELP         GZMM 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        GPC3         TFRC         EGR4      LRRN4CL         OSR1     MIR155HG        THBS4      GADD45G     APOBEC3G          OSM        SOCS1     C16orf89      TSPAN13 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
     TM4SF18          EZR       TUBB2B        OLFM1      ADAMTS5        LAMP5     HIST1H1D        ITM2C          TRH    TNFRSF13C     MAP3K7CL    LINC02446     TMEM176A 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        NRGN        FIBIN         ENHO         NCR3         DGKI     MARCKSL1     TNFRSF17        RSPO3         BDNF       COL4A2   AC104024.1      TMEM215        MYOZ2 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
        CD3E           F5       CHMP1B       IGFBP1     IGHV3.53       CD40LG     IGKV2.30        CFHR1        HBEGF         LAG3        DUSP4         BATF 
0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 
 [ reached getOption("max.print") -- omitted 2000 entries ]

# To prepare the long-format data:
shap_long <- shap.prep(shap_contrib = shap_contrib, X_train = ds2_train_data$data, top_n = 20)


# is the same as: using given shap_contrib
# shap_long <- shap.prep(shap_contrib = shap_values$shap_score, X_train = ds2_train_data$data)
# (Notice that there will be a data.table warning from `melt.data.table` due to `dayint` coerced from
# integer to double)


# sometimes for a preview, you want to plot less data to make it faster using `dilute`
shap.plot.summary(shap_long, x_bound  = 1.2, dilute = 30)

## FB
shap_contrib <- predict(reduced_ds2model, ds2_train_data$data, predcontrib = TRUE)
shap_contrib <- data.frame(shap_contrib[[4]]) ##提取SMC1对应的分类
BIAS0 <- shap_contrib[, ncol(shap_contrib)][1]
shap_contrib[,"BIAS"] <- NULL
imp <- colMeans(abs(shap_contrib))
mean_shap_score <- imp[order(imp, decreasing = T)]

shap_values<- list(shap_score = shap_contrib, mean_shap_score = mean_shap_score, 
    BIAS0 = BIAS0)

# To prepare the long-format data:
shap_long <- shap.prep(shap_contrib = shap_contrib, X_train = ds2_train_data$data,top_n = 30)

shap.plot.summary(shap_long, x_bound  = 1.2, dilute = 30)

shap.plot.force_plot(plot_data,zoom_in_location = 800)

Data has N = 6675 | zoom in length is 150 at location 800.

shap.plot.dependence(data_long = shap_long, x= "LUM",
                     y = "BGN", color_feature = "BGN")   

`geom_smooth()` using formula 'y ~ x'

shap_int <- predict(reduced_ds2model, ds2_train_data$data, predinteraction = TRUE)

Error in predict.xgb.Booster(reduced_ds2model, ds2_train_data$data, predinteraction = TRUE) : 
  std::bad_alloc

data("iris")
X1 = as.matrix(iris[,-5])
mod1 = xgboost::xgboost(
  data = X1, label = iris$Species, gamma = 0, eta = 1,
  lambda = 0, nrounds = 1, verbose = FALSE)

# shap.values(model, X_dataset) returns the SHAP
# data matrix and ranked features by mean|SHAP|
shap_values <- shap.values(xgb_model = mod1, X_train = X1)
shap_values$mean_shap_score

Petal.Length  Petal.Width Sepal.Length  Sepal.Width 
  0.62935975   0.21664035   0.02910357   0.00000000 

shap_values_iris <- shap_values$shap_score

confidence曲线

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.